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ABSTRACT 


The solution of the nonlinear two-dimensional reactor dynamics 
equation subjected to prompt feedback conditions using the finite 


element technique leads to the matrix formulation AsV, - See on 


Coals (i,j,k = 1, ..., N). This system has been solved directly 


in a previous work; but because the nonlinearity C was 


ijk" 5 "k 
premultiplied by Ee large computer storage was required for the 
small problem considered. The task of this thesis is the development 
of computational techniques which allow the problem to be solved for 
large systems. Specifically, these techniques are: (1) the treatment 
of the nonlinearity on the element level, (2) the compacting of the 
sparce matrices to include only non-zero terms, and (3) the construction 
of a new computer code based on the Crank-Nicolson formulation for the 
solution of differential equations. 

To support tne theory presented, test problems were solved by the 
original method, the linearized technique, and the Crank-Nicolson 
treatment. The results were analyzed and compared graphically. All 


three of the innovations developed in this thesis appear to be useful 


tools for solving nonlinear time dependent differential equations. 








TABLE OF CONTENTS 


ie INTRODUCTION -------------------------------------~------------ 9 
II. REDUCTION OF THE NONLINEARITY --------------------------------- 1 
A. GLOBAL TREATMENT OF THE NONLINEARITY --------------------- 14 
B. TREATMENT ON THE ELEMENT LEVEL -------------------------~- Va 
III. SOLUTION PROCEDURE WITH THE LINEARIZED FORM ------------------- 19 
A GENERAL TREATMENT OF THE MATRIX EQUATION ----------------- 19 


B. PRACTICAL CONSIDERATIONS WITH THE LINEARIZED TECHNIQUE --- 20 


C. DIFFICULTIES WITH THE APPROXIMATION ---------------------- 21 
IV. THE Nxp COMPACTING SCHEME ------------------------------------- oP. 
A. ESTABLISHED COMPACTING METHODS --------------------------- 22 
B. BASIS FOR THE COMPACTING SCHEME -------------------------- 23 
C. CONSTRUCTION OF THE Nxp ARRAY ---~------------------------ 24 
v. THE CRANK-NICOLSON FORMULATION ---------------- sees 26 
se Ore PAR EPeMern@ns =—-———————————---—---_-_-__________- 30 
eee ie (eneee hhc AND RESULTS =-——————————===—=——— === 31 
A. THE PHYSICAL MODEL --—-------------------—----—--_____-___.. 31 
B. COMPUTER PROCESSING CONSIDERATIONS ---------------------~- at 
G5 PRORD EM AOU C 1S SS ee 31 
D. EVALUATION OF ERROR ---------------------------+--~------- 34 
VIII. CONCLUSION ---------------------------------------------=----=- 35 
APPENDIX A — TABLES ~--------=—=—-----———-----—=~——-—-------—------—- 36 
CUP. ATCO ee 39 
APPENDIX C - COMPUTER PROGRAMMING CODES ----------------------------- 61 
LIST OF REFERENCES --------------~---------------~--------—----------- 120 
INITIAL DISTRIBUTION LIST -----~-------------- ae SS 121 








10. 


lee 


Tie 


Juste 


14. 


key 


lular 


ET. 


18. 


oe 


20. 


Za. 


List OF FIGURES 


REACTOR MODEL WITH 38 NODES ------------------- 


NODAL CONNECTIVITY ---------------------------- 


REACTORSNMODEL Wide 132 NODES ------------------ 


PLOGsObeCEmeEe Dp ESEURBANCE = CENTER POINT ----- 


COMPARISON - 


CENTER DISTURBANCE —- CENTER POINT 


PEGieor eCENBER DISTURBANCE — CORE POINT ------- 


COMPARISON - 


GENTER DESTURBANCE - CORE POENT -= 


POT OF CENDER DISTURBANCE = REFLECTOR POINT -- 


ae ee ee eS SS ee ee ee oe 


= SS es ee ee eS ee Se 


COMPARISON - CENTER DISTURBANCE - REFLECTOR POINT ------------- 


PLOT OP avi ORM DES TURBANCE = CENTER POINT ---- 


COMPARISON - UNIFORM DISTURBANCE - CENTER POINT 


PLOT OF UNIFORM DISTURBANCE -— CORE POINT 


COMPARISON - 


PLOT OF UNIFORM DISTURBANCE - REFLECTOR POINT 


UNTFORM DISTURBANCE - CORE POINT - 


ea eee ee i 


SSS ee ee ee ee ee oe ee 


COMPARISON - UNIFORM DISTURBANCE - REFLECTOR POINT ------------ 


PLOT OF SKEW 


COMPARISON - 


PLOT OF SKEW 


COMPARISON = 


PLOT OF SKEW 


COMPARISON - 


DISTURBANCE - CENTER POINT ------- 


SC EWeDEStUREANCE: = "CENTER POINT -- 


DISTURBANCE - CORE POINT -------~---------3-------- 


SKEW DISTURBANCE - CORE POINT ---- 


Dis tRBeaNCE — REELECTOR POINT ---- 


SKEW DISTURBANCE - REFLECTOR POINT 


— a ee ee ee Se ee ee ee ee 


40 
41 
42 
43 


44 


46 
47 
48 
49 
50 
Syl 
2 
53 


54 


56 
a 
58 
op, 


60 





ACKNOWLEDGEMENTS 


With grateful appreciation I wish to acknowledge the persistent 
and dedicated efforts of my thesis advisors, Associate Professors D. 
H. Nguyen and D. Salinas. 

I sincerely thank my many friends on the C.W. Church Computer 
Porgramming staff who contributed significant advice and direction 
toward the successful construction of my projects. Paramount in this 
group was Mr. Roger Hilleary whose mathematical expertise coupled with 
his detailed knowledge of available software clarified my thinking and 
provided the necessary insight for the task at hand. 

Frequently there exists a noticable gap between the knowing how 
and the getting it done. The people who were able to help me close 
this gap were the Computer Center Operators, a really fine group of 
professionals. Among them was Mr. Ed Donnellan, a very patient man 
who led me through many intricacies of the job control language. To 
these individuals I owe much gratitude. 

I appreciate my opportunity to study for the past two-and-a-half 
years at the Naval Postgraduate School. I thank those professors with 
whom I have worked for their contribution to my total educational 
experience. This experience includes not only an increased technical 
expertise but also a deeper maturity and a broader horizon. Among 
others, the names of T. Sarpkaya, J. Brock and T. Cooper will long 


stand out in my memory. 








Lt. INTRODUCTION 


When temperature-dependent feedback is considered in the nuclear 
reactor dynamics problem, a nonlinear field equation in space and time 
results [Ref. 1]. For the non-homogeneous, or multi-region reactor, 
however, the space-dependent dynamics behavior following a nonlinear 
initial disturbance is no longer reachable in analytical form. 
Fortunately, by modeling the reactor as a system of finite regions of 
interest where the neutronic properties are known, the method of finite 
elements can be applied to yield the solutions. The fundamental concepts 
relating reactor behavior to the finite element formulation have been 
presented in 1974 [{Ref. 2], but a more recent work by Nguyen and 
Salinas [{Ref. 3] gives a thorough discussion of the complete problem. 
That work forms the basis and starting point for this thesis, the 
objective of which is to develop improved computational methods for 
dealing with that finite element formulation. 


In Reference 3 the dynamic flux equation under prompt feedback 
conditions was given as 


ap (r,z,t) 
m = Ze - 2 2 ; 
dt i ve Po Vea a ho oY (1) 


a 
V mam'm mm am’m 


for each non-homogeneous zone "m" contained in the reactor body, 


where the usual symbols are employed: 


= axl 
han Vv ee a 


K =e Bian! Ose (°C/unit flux-sec) 


-] 
— f 
a (A/V er Ce) sec 





p a = heat capacity 
a@ = reactivity temperature coefficient ECR 
h = convection heat transfer coefficient 


A/V = heat transfer area to volume of energy 
generation ratio 


de = neutron fission cross-section 
% = neutron absorption cross-section 
vy = neutrons emitted per fission 


e = energy produced per fission 
The subscript 'm' will be omitted from further discussion for simplicity. 


The transformation of this general equation into a problem in 
finite dimensional vector space begins by constructing the following N 


term approximation 


bad 


N 
VeGE gic =) 7) (EGG) (2) 
Ld. ] 
J 
where N is the number of degrees of freedom, or nodes, in the space, 
and the oe are the basis functions of the approximate solution space, 


The Galerkin method seeks to make the residual 


~- 


R(x,t) =Ly- f, 


where Ly = f is the field equation, othogonal to each of the basis 


functions: so that 


1 ; : 2 , ; : 

This discussion leading to the establishment of the matrix 
equations is essentially an abridgement of the development given in 
reference 3. 








[S,@R@tde=0 151, 2..5N, (3) 
S 1 


For the nuclear reactor dynamics problem of Eq. (1), the residual 


becomes 
Ris.t) = ay" —- V bys ~ VAE bp + wi fe (4) 
—s dt a a 
with VoRS =w . Putting this expression into Eq. (3) and integrating 


by parts (a distinct advantage of Galerkin), the following coefficients 


emerge 


A j “Sf ex rdrdz (5a) 
S 








ve Be WG. BE 
= a , Jk A 
oe we ancy mg ee) 
Cray = SL GG 5G, rdrdz (5c) 
i Si a 


in cylindrical coordinates where dS = rdrdz. This allows Eq. (1) to be 
written as follows for the case of uniform neutronic properties within each 


reactor region Mm. 


N 
3 AL Ys =- VD) » Te ee roel >, Ay s¥5 
(6) 
N 
J=1 K=1 
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This becomes, upon combining constants, in Einstein summation convention, 


where there is no sum on an underlined repeated index, 


A. b. = AB -~ w.C 


Eins fio ye rene gree 


= (7a) 

¥, = £,(,,¢) ig mos ey i 
With boundary conditions, the system is now well posed and ready for 
solution. The equation-solver must be chosen with care, however, 
because this system of ordinary differential equations is both stiff 
and non-linear. "Stiff" is used here and throughout this thesis to 
denote a system which gives a large response to a small stimulus; in 
this case, the nature of the reactor dynamics problem predicts a large 
change in flux over a small change in gee? Current experience 
suggests Gear's predictor-corrector method [Ref. 4] as written in the 
computer programming code DVOGER [Ref. 5]. It requires only 1) the 
locus (value) of the points at a given time, 2) a routine to evaluate 
the instantaneous derivatives at any time and 3) a routine to evaluate 
the Jacobian (J = 3£ Cyst) /3y,). Based on this information, the program 
assumes a time step, which may be quite small (the minimum size being a 
munetion Of the particular computer), and tries to fit a trial (predicted) 
solution for that time interval. The corrected (integrated) solution is 
then obtained by iteration until convergence is attained. If the 
predicted and corrected solutions fail to match within a specified error 
criterion, the time step is decreased and the process repeated. On the 
other hand, if "excessive" accuracy is attained, the next attempted 
time step will be increased. Gear's method in DVOGER is well suited 


to stiff non-linear systems. Aliso contained in DVOGER is the Adam's 


Paes ae 1 
This definition is somewhat broader in scope than that used by 
IMany other authors. 


ina 





method which parallels that of Gear except that the Jacobian is not 
computed. Lacking the additional information about the rate of change, 
the Adam's method must proceed much more cautiously with stiff systems 
and hence marches forward with exceedingly small time steps. Present 
experience confirms the caveat advice given by DVOGER that the Adam's 
method is not intended for "stiff" cases. 

The difficulty that arises from acquiring solutions to Eq. (7), 
then, should not, and, in light of the present experience, does not 
result from DVOGER or any other acceptable equation solver, but is, rather, 
a function of matrix size, vis a vis the number of nodal points in the 
finite element approximation. The finite element modeling of large 
reactors or the close scrutiny of small ones is severely restricted, 
therefore, unless the number of nodal points can be raised significantly. 
As an example, the test case considered in this thesis consisted of only 
38 nodes but, when processed directly from Eq. (7), required over 300K 
bytes of storage in the IBM 360/67 computer. This requirement came 


largely from the ra) ttc ] matrix alone needing a size of 4(38x38x38) = 


ok 
219K bytes in single precision. Manipulations with this cubic were 
indeed quite burdensome. Since machine processing time is, other things 
being equal, dependent on size, the combined time and space requirements 
would prohibit the consideration of even moderate sized problems except 
on the largest computers. Additionally, the inversion of [A] in Eq. 
(7) is indicated in order to provide an explicit value of " for the 
DVOGER routine. It is this difficulty of size, and with it processing 


time, to which the remainder of this thesis is devoted. Experience with 


the DVOGER routine is presented, and the authors equation-solver based on 


ite2 





a Crank-Nicolson formulation is discussed as the analysis of the size 
problem is developed. Table I shows clearly the influence of matrix 


size on computer storage requirements. 
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II. REDUCTION OF THE NONLINEARITY 


A. GLOBAL TREATMENT OF THE NONLINEARITY 
The evaluation of Eq. (7) generates the following global matrix 
system, where the subscripts indicate the matrix size required for N 


number of nodal points 


| | {Yara f : ea {xsi} ae esx Yr) 2} (8) 


where w is the time response of the flux at each nodal point and re 

is the product Le (jai ae Ne coolved in this form, the major 
limitation of the procedure was the large ne size of the cubic array 
(ewiten is introduced as the result of the nonlinearity of Eq. (1). 

With this size requirement rapidly increasing for a finer mesh consisting 
of more nodal poring the investigation of Reference 3 was limited to a 


small but practical N of 38. 


B. TREATMENT ON THE ELEMENT LEVEL 

The cubic form of C need not appear if it can be shown that the 
nonlinear portion of Eq. (8) can be broken apart into the product of 
two linear components such that one component can be evaluated on the 


element level and the element contributions then summed into global form, 


In general, for a nonlinearity of order m, a finite element 
problem of N nodes will generate a matirx size of ym, 


14 





resulting in a “regular' NxN matrix C* ready for multiplication by 


the remaining w also summed to the global level. Proof that this is 


indeed the case, i.e., 


| Sx nant ¥en 2 } - [a | {Yn } (9) 


is outlined as follows: 


Reference 3 established the element nonlinear term as 


Sf Og54) OE IAA = Chay VACE) H(t) (10) 


= 


where the integral is over the element area A. ; 


eS 


oi Ah, G _ardz, ieee 2.5 (11) 


and ¢ being the local triangular coordinates, an innovation of Felippa 


[Ref. 6]. The coefficients 
a regularity which suggests 


of this integral results in 


Alstah 


(e 
Il 


iy 


o) 
I 


2) 
ll 


Cc) 
tl 


QO 
ll 


O 
Il 


‘o) 
Il 


a ai form a 3x3x3 element array possessing 
a kind of "cubic symmetry.'' The evaluation 
the following expression 
y°[24r, 4- 6(r, + r4)| 
y°[6r, + 4X, + 3x7 ] 
y°| 6x, + 2x, 3F Ars | 
y°[4r, + 6r, + ar, | 
7 ae + ty +- r)| 


Y “Lar, + 2X, + 6r , | 


crZ) 


Y “lor, + 24r., 3 6x. 


ey + or, + 24r | 


Ibs: 





where 
y = TA ,/180 


and 


Sig 7 ai oes 


7 es ee OS ate asidneRef. 3. 


Assuming that Ue may be broken apart as Wes = ww * ¢, the LHS 


(left hand side) of Eq. (10) is rewritten as 


Lf i G¥p Ged 4,5,K = 1,253 (13) 
i 
e 
or 
y b 
21 Lf [615553] | ¥ | [6,555] | o, | © dr dz (14) 
A ¥3 %3 


e 


Expanding and collecting terms yields 


3 2 2 
ai // oe |, Joy ey5) + robyoo + 630453) 


e 
se hs (ie 2 + 2 + ) 15) 
ee pt “seo ( 
+ ¢.(r 5 + Y + x a 
a ees eco 


2 2 2 3 2 
( Yr 
Foy $4845 ySo + 69550 + 63by5053) + $50r] 5450 + Toby + 145054) 
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Z 


2 
+ $307) 0)5554 + To554 + 130555) 


2 2 2 2 
PVs 14 %y 5453 F Fobq5b3 + 655,53) + b9(r 1245005 + Tobob, t 140505) 


2 2 
we gS Ey oven) 39} | dr dz a ra! 


Since the quantities v. and . (i,j = 1,2,3) depend only on time, 
they are unaffected by the integration. Further, the expressions 
inside the parentheses are given a unique name are EO fixeetiews 
position in relation to a particular ¥595 for some k= 1,2,3. Now 


the above expression may be written as 


21 ¥59, If disk dr dz , (16) 
A 


Using the integration formula from Felippa [Ref. 6] 
Q m on = g!'m!in! 
Wh Sty ea) ee Gey * oe C7) 
A 
e 


where &,m,n are the exponents of ¢f in a specific d. The integration 


5S) 
of (16) gives 


(1A /180) Vio ere (18) 


and f.. 
a 


a - : 
jk * TA ,/180 is identically equal to Ci sk Of (12) 


The indicated summation is conveniently expressed by 


De fit =o os See 2 (19) 
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which is immediately recognized as moves the desired 3x3 


e 

hk 
k 

element matrix. In this form it is combined into the global [C.. oy! » 


and the linearization has been accomplished. As will be shown later, 


it may be convenient to let 96 = jp when for small At, $6 = 


toc 


aa ~ De . This approximation is useful for predictor-corrector 


equation solvers in order to allow construction of [C only once 


* 
Nx! 


for each successful time step. 
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III. SOLUTION PROCEDURE WITH THE LINEARIZED FORM 


For clarity, throughout the remainder of this thesis the following 
naming convention will be used: x? denotes the nonlinearized direct 
treatment by DVOGER as solved in Ref. 3. "Linearized approximate" 
indicates the linearized technique based on element level multiplication 


Or C by w The "linearized exact" method, however, is the element 


t-At- 
level multiplication of C by the trial ve The compact, or Nxp 


formulation may be used with any of the linearized methods in addition 


to the CRANKO treatment discussed later. 


A. GENERAL TREATMENT OF THE MATRIX EQUATION 

Solution of the matrix formulation of Eq. (8) involves clearing 
the LHS to provide a discrete relation for each b- This is accomplished 
either by iteration, which yields an approximation, or by matrix inversion, 
which has the advantage of giving the individual by explicitly. In 
both the nonlinearized nw? and the linearized NxN forms, Eq. (8) is 


eee eater itcation byw Ts) =; thus 

fa] TLA](b} = [AD “[B]{y} + [al “Ece){y} (20) 
in the linearized ye = oo) case or as 

[AD {AN (G} = (A) “ EB] {y} + CA) TIC) {y*} (21) 


tor the direct, nonlinear, formulation. This results in 


- 


by SMGNSILG su Pee) Con (22) 





where 
1@ 1-2 (eS on. (Kel (23) 


as appropriate. The point here is that the [AB] and the sy are 
formed once and for all based on the time independent properties and 
geometry of the problem, whereas the [Cd] must be computed for each time 
t. Although the construction of the [C] is required only once in the 


direct (n>) method, it must be additionally multiplied by ee alte 


each time step. 


B. PRACTICAL CONSIDERATIONS WITH THE LINEARIZED TECHNIQUE 

Relating Eq. (22) to the demands of the subroutine DVOGER, which 
involves the evaluation of w for every trial value of yp, it is 
apparent that the construction of Tene ccna’: fOmreachsot thee vem s 
could be a time consuming process, and hence a serious drawback of the 
linearized treatment. Noting, however, that for small At, yp = ane 


the approximation [Cv] = [Cv], where $¢= yp » overcomes this 


C= 5c 
disadvantage by allowing the construction of iP iat only once for a 
given time t, regardless of how many trial solutions are attempted by 


DVOGER. Equation (22) can now be further simplified since [A] + Ico] 


is readily combined with [AB] to reformulate the problem as 
: kk 
{vw} = [Cc Jf{y} . (24) 


Thus, although the linearized form can be handled exactly, i.e., with 

: ke 
¢ = ~, it may be more expeditous to establish the [C ] only once per 
valid time point, (that is, where the convergence criterion of the equation 
solution has been satisfied). In this case, the linearized treatment is 


then an approximate technique in that 9 = eee # ve 
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C. DIFFICULTIES WITH THE APPROXIMATION 

The simplification expressed as Eq. (24) works extremely well when 
used with DVOGER as long as large changes in wy are incurred by small 
At (the so-called "stiff" system). When the solution approaches steady 


state, however, approaches zero faster than (yj )/At does, 


t  Yt-at 
and this causes a second kind of perturbation wherein the system becomes 


increasingly sensitive to small changes in y Apparently 


t  Yteat’ 
DVOGER requires very accurate " information. The W(t) provided from 
Y(t) x CLONE) does not match the DVOGER prediction of Yeey Which is 
based on the assumption that Yee) was generated by Wee" This 
discrepancy in the neighborhood of 7 =Q causes DVOGER to be unable to 
recognize the steady state solution. Prepared especially for stiff 
systems, the equation-solver expects, instead, a large scale change in 

y. The time step is therefore reduced accordingly, which means that very 
slow progress is made in the steady-state region which is otherwise 
handled quite rapidly under the exact methods. The implication here is 
that, in those cases which reach terminal flux values rather early in 
problem time, the exact formulation of the linearized treatment is to be 
preferred, even though the rc must be computed for each attempted 
ean This aspect has been analyzed for the test problems considered, 


with the results given in Table II. 
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IV. THE Nxp COMPACTING SCHEME 


A. ESTABLISHED COMPACTING METHOD 


The problem under consideration here, being of the form Ay iY, 
Bay ree Vy)» is formally handled by multiplying through by lan 


discussed earlier. For large matrices, the inversion process, and indeed 
all matrix operations, become extremely costly in terms of machine time 
and storage. The fact that the matrices may be banded or symmetric does 
offer significant economy when the matrix is not inverted. For example, 
a symmetric NxN matrix may be stored as NxN/2. For matrices resulting 
from finite element formulation, the numbering schedule determines a band 
width to be used with conventional compacting schemes. In the finite 
element case, the difference is compared for each node in the mesh. The 
largest of these differences plus one equals the minimum band width 
allowed. For example, node 7 in Figure 1 has a maximum difference of 

13 — 2 = ll, whereas node 22 has the largest difference (37 - 10 = 27) 
for the system. Thus the resultant minimum band width would be 28. As 
the mesh increases in nodal points, the minimum band width must also 
increase in response, no matter how adroitly the numbers are assigned. 
Regardless of how sparce a matrix may be, these schemes all ensure that 
the ''compactness'' grows commensurately. In addition, the calculation of 
an inverse becomes ever more unpleasant, especially since this operation 


on a sparce matrix generally results in a dense, nonsymmetric inverse. 
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Be BASIS FOR THE COMPACTING SCHEME 
In solving a matrix equation Ai 4%; = Be (y,> ake ae vy eats Tot 
necessary to form the inverse, but rather, multiply out the left hand 


side and rearrange to obtain a set of equations: 


where i,j,k =1, ...N. This procedure is especially suitable for 
sparce systems since the number of non-zero ay entities will be small. 
Further, if these non-zero oe can be located and tagged with an 
identifier, then they may be stored together in a dense array, thereby 
replacing the standard NxN size matrix by an equivalent array of size 
Nxp, where p is the maximum number of non-zero entries on any one row. 
This type of compacting is especially well suited to matrices 
resulting from finite element formulation, with p determined by the 
choice of the finite element and the discretized model. In the case 
under study, a linear triangular element is used such that, when assembled 
on the global level, it forms, for the purpose of illustration, a network 
of interlocking hexagonal polygons, each having an apex located directly 
over its parent central node. Such a pattern is sketched as Figure 2. 
Note that the sides are indeed the contributing "neighbor" nodes. Boundary 
nodes will, of course, be missing any "would be" contributors sought in 
the boundary exterior. In a regularly drawn system, as in Figure 2, an 
interior node will have six neighbors plus itself for a total of seven 
contributors. That is to say, that no matter how large the system of 
mesh points, a relation describing the effect of the system on any given 


point would contain a maximum of seven non-zero values. 
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Geometric considerations or a desire to examine some area of the 
mesh more closely may result in systems of irregular discretization 
wherein the maximum number of entries may be greater than seven. This 
analysis may be extended to consider other finite element shapes in a 


a?milar manner. 


C. CONSTRUCTION OF THE Nxp ARRAY 

To utilize this attribute of finite elements to produce the dense 
Nxp system, it is necessary only to construct a table of nodal points 
listing their contributors. This table, or connectivity array, is then 
used as an index to locate the values of the various quantities associa- 
ted with a particular node. All the standard matrix operations can be 
performed on these compacted arrays, but the construction of a matrix 
inverse has no usefulness since a NxN system is then reconstituted from 
an Nxp array. 

As an example, construction of the connectivity vector associated 


with node 18 in Figure 1 is 


felicer) 27 eon 9 4 13 12) 


and the vector associated with node 6 is 


ioe wel 6. 0-0 6). 


The contributors may be entered in any order except that, for computa- 
tional ease, the central, or "parent" node is the first element. These 
vectors are collected into the node neighbor connectivity array of size 


inp, Or, for Figure 1, 38xs. 
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To multiply a Nxp = array of as; elements by a vector of 
elements, a search of the connectivity array is performed on the "i"th 
row to locate a match for the Ber When the match is found, the values 
represented by a, and oe are multiplied together resulting in a new 
C, vector. If the system is correctly formed and compacted, there 
always will be a match. 

The compacting scheme is bound by the geometry of the problem, and 
as such, has meaning only as it pertains to that problem. For example, 
operations with two equi-sized Nxp arrays of different connectivity 
relations cannot be perfermed. This is in contrast to regular matrix 
Operations where the restriction is only to size and not to origin. A 
short computer program to demonstrate the operation of this Nxp-_ scheme 


is given, with results, in Appendix C. 
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V. THE CRANK-NICOLSON FORMULATION 


Recognizing the very sparce nature of the system of differential 
equations formulated by this problem and, if the matrix inversion process 
is to be avoided, the necessity of iteration, a straightforward attack 
based on the definition of the derivative (finite difference) appears as 
a feasible alternative to the DVOGER technique. One such approach was 
presented by J. Crank and P. Nicolson in 1947 [Ref. 7] and has been 
discussed in many works throughout the ensuing years. This method has 
been shown to be unconditionally stable. An adaptation of this method 
is given here as follows: 

In the general sense, a matrix formulation of a set of linear 


differential equations, can be represented as 


= + is] = eee 


where w is a function of time. For the case of an initial disturbance 
emeinestorcing fLumetion, F(t) = 0 for t > 0. Focusing attention for 
the moment on a particular equation and expressing [16] in terms of the 


fetinnelon Of the derivative 


py p ty 
Eom AB || ot t b28e 
A are || Aas C Po (27) 
wheaeh is 
fA = 1/2*At-C]p, = {Areeel/ 2-AceGy yy 4s (28) 
[(2/dt+A)-C] yp, = [(2/At-A)t+C])_ A, | (29) 
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letting D = (2/At°A)-C, E = (2/At*A)+C, 96 = Vea then, returning 


te 


to the system form the result is 


D..wW, = Bo oe = E, 53) = ibe eeters LN Se (30) 


This constitutes a set of linear simultaneous equations which may be 
solved by any convenient method. Chosen here is Gauss-Seidel iteration 
[Ref. 8] which is well suited for the problem at. hand since, the sum of 
hay and Baar will contain very few non-zero terms, thereby incurring 
relatively small roundoff errors regardless of system size. 

This procedure, which also uses the Nxp compacting, has been 
written for this thesis as the computer code CRANKO. After receiving 
the matrix information and control parameters from the calling program, 
CRANKO first builds the C matrix cf Eq. (26) for the initial ¢ value. 
Then it selects a trial time interval and, forming the equation set (30) 
based on this At, attempts to find a solution satisfying a specified 
convergence criterion based on a relative error. If no satisfactory 
solution set is found after a number of iterations, a smaller At is 
selected; and, after recomputation of Eq. (30), another attempt at 
convergence is made. The successful solutions then replace the previous 
y thereby forming the starting point for another cycle, beginning with 
anew computation of the C matrix. The choice of At is the control- 
ling factor. The scheme employed counts the number of iterations required 
for a solution. When an increasing number of iterations indicate greater 
difficulty, the interval At is decreased; whereas, if fewer iveraene 
are required, the time step may be relaxed and made larger. Experience 
with the CRANKO routine establishes this method as an extremely fast 


technique for solving this type of problem. 





For those readers unfamiliar with the Gauss-Seidel iteration scheme, 


a simple example using a 3x3 system is given: 


di, 4y2 gg al 
doy don gg 1 | My | = | 22 
dey dan dag | LY3 = 


me-ct tteration: 
dyy¥y + dines F 41 3%5 = &y 
¥, = Ley - Cdyody + 4) 36,)1/4), 
* 
YF YY 


* 
Bag) Seg hiey Cee = er 


* 
vy = ley — (dyi¥y + 45394) ]/ dy, 


x 


Ye 
| it eae eae E 
31¥1 + F30¥2 t 453%3 = &3 
= fae a an yaa 
v3 = fe, - (daiv) + 4,505) 1/d,, 


for successive iterations, 
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* * 
fag gee NA ay 


2 ee ed. eye 
1 = fey ~ (dy ov. + dy 3%) 1/4), 
KK 
veal 
eK * _ 
doi¥y + Fogo + 44303 = & 
K* * 
Vo = [ey ~ (doi¥y + 4,393) )/4,, 


prx_ vy 


27 2 
q eck q Od F _ 
31% + 430¥2 F 433%3 = &3 


= 7 K* 
v3 = [e, - (dg¥) + d35¥. D1/d3, 


etc. 
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VI. SUMMARY OF AVAILABLE METHODS 


For clarity, the innovative computational tools introduced in this 
thesis are summarized. First is the linearizing of the governing equation 
Eq. (8) by the multiplication on the element level of Cr ak by De 
From this development, two options appear feasible: either the exact 
computation of [Cc] for each trial yw, or the approximation of 
ve SE to allow the construction of [c"] only once during the 
search for a new v.- This approach uses, as does the "direct," or 
menlinear treatment involving [C], the inverting of [A] to clear the 
mes as shown in Eq. (22). 

The second innovation is the reduction of the NxN matrix to a 
smaller Nxp compact form based, not upon the mathematics of the problen, 
but upon the geometry of the finite element model. The same solution 
methods may be used, with the only difference being that the inverting 
of [A] is replaced by iteration. This iteration can only be as accurate 
as the imposed error criterion. For large systems, however, where the 
Al im produces a full NxN matrix on the RHS, this iteration is 
especially beneficial. 

The last contribution was the development of a new equation-solver 


based on a different theory than the DVOGER routine. The new code deals 


with the differential equations directly and employs the Nxp compacting. 
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VEL. TST PROBLEMS AND@RESULTS 


A. THE PHYSICAL MODEL 

A cylindrical nuclear reactor consisting of a core region and 
blanket having overall a height of 220 cm and a diameter of 180 cm was 
modeled for the test problems. As shown in Figure l, a radial slice 
from this reactor was discretized by a finite element of 54 elements and 
38 nodal points. After a detailed analysis of this mesh was completed, 
preliminary investigation was begun using a model of 220 elements 


(Figure 3). 


B. COMPUTER PROCESSING CONSIDERATIONS 

All programs were written in the FORTRAN IV language and processed 
on the IBM 360/67 computer using the FORTRAN 'H' compiler. The results 
presented in terms of CPU storage requirements and processing time ought 
to be indicative, in a relative way, of those obtained from other machines. 
Single precision (seven significant digits) was used throughout the 
investigation, with the results later verified using double precision 
(fifteen significant digits). The observed relative discrepancy was less 
than .01%. It is recongized, however, that for larger systems a more 


significant difference may result. 


C. PROBLEM ANALYSIS 

The three techniques discussed have been incorporated into computer 
codes and have been tested on those problems described in Reference 3, 
namely, the dynamic reactor response to 1) a disturbance at the center, 
2) 2 uniform disturbance throughout the core, and 3) a disturbance occurring 


3 
at a skew point in the core (R= 40 cm, Z= 0m). The direct, or N 


oH 





Matrix (an Nx (N2) array) equation was solved in Reference 3, and those 
results are used in this thesis as a standard for comparison. The 
computer codes developed here have been constructed within the framework 
and parameters of the original nN? program, thereby allowing any 
discrepancy to be attributed to the technique employed. For each tech- 
nique and disturbance, the flux at three sample points was recorded 
throughout the time history of the solutions and written onto a computer 
storage device. A separate program was written to process this infor- 
mation and produce the correlations in tabular and graphic form. The 
graphic results are presented in Figures 4 - 21 for the center point 
(k= QO9em, Z = 0 cm), a core point (R = 40 cm, Z = 40 cm), and a reflector 
noimi (kh = 7/> cm, Z = 80 cm). Table IL lists the specific programs 
employed arranged by disturbance type and gives parameters such as 
storage requirements and computer processing time. 
ANP Solutions with the CRANKO Subroutine 

All techniques appear to give acceptable accuracy, however 
there is a startling difference in processing time. The CRANKO subroutine 
was able to accomplish in seconds the work which required minutes for the 
DVOGER subroutine (See Table II). The steady state solution given by 
CRANKO was in every case higher than that obtained by DVOGER. The 
implication here is not clear since there is no analytical solution 
available to represent the exact answer. In general, the transient 
solutions agreed well, except for the oscillation of the CRANKO solution 
in the early stages, most clearly depicted in the extreme case of Figures 
8 and 9. This fluctuation is most likely the result of inadequate control 
of the CRANKO time step selection where apparently an interval has been 


chosen which is too large to maintain good accuracy. It is noteworthy, 
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however, that the technique is so stable that even after large deviation 
it is able to recognize the error, correct itself, and return to the 
solution curve given by the DVOGER routine. 
2. Solutions with the DVOGER Subroutine 
The linearized technique using premultiplication by the previous 


time step flux (v ))contormed to tWeoretical expectations. As the 


Eee 
steady state solution was approached, this method proceeded with in- 
creasing difficulty. Hence much more processing time to advance through 
the steady state was required. 

The direct ne method, on the other hand, has no difficulty 
in recognizing the steady state and marches forward with large time 
increments until problem termination. A very successful effort was made 
to use the exact treatment with the linearized technique. As expected, 
more processing time was consumed in getting to the steady state due to 
the requirement to recompute the [Cv] for each trial wy, but once the 
steady state was approached, the procedure progressed rapidly as in the 
direct treatment. Comparing the direct with the linearized exact method, 
theory developed in this thesis predicts identical results with twice 
the computing time required by the linearization. As shown in Table II, 
this is indeed the case. The largest relative error observed between the 
two methods was .06%, which is attributed to roundoff steming from the 
doubling of the computational requirements. When the greatest accuracy 
is required, the exact treatment with linearization may therefore prove 
the most useful, especially if the problem rapidly reaches steady state. 
This discussion of the linearized technique is applicable both to the NxN 
method or to the Nxp compact method. However since the Nxp method 
requires iteration, more processing time and less accuracy might be 


expected. 
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D. EVALUATION OF ERROR 

The error analysis for these procedures must remain subjective as 
there is no "correct" result. Several sources of error are present and 
deserve comment. Tests were made with both EPS = .1 and EPS = .01, the 
error criterion for the DVOGER routine, and little difference was noted 


for the large increase in processing time. [EPS is defined by DVOGER as 
; ERROR (i) 
i=l “MAX 


The convergence criterion [relative deviation of any 


v. between successive iterations] for the iterative solutions in CRANKO 
was varied by tenths from .01 to .001 with the result of producing higher 
steady state values of only a few percent as the criterion was made more 
stringent. Processing time for this variation increased only a few 
seconds. Double percision trials for the linearized DVOGER method and 
the compact CRANKO were made with neglibible difference. 

The grid of only 54 elements theoretically introduces the largest 
source of error from the "true" solution. With the techniques developed 
in this thesis, a finer grid size of 220 elements has been implemented 
into the CRANKO and the linearized codes. Because the linearized code 
(DVOGER) requires about 566K core size, extensive testing with this 
technique is not anticipated, but the CRANKO routine which requires only 
188K is currently being investigated. As the results of this investiga- 
tion become clear, a better understanding of the error induced by each of 


the various approaches to this problem would be obtained. 
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VIII. CONCLUSION 


It appears that nonlinear differential equations resulting from a 
finite element formulation (or for that matter finite differencing) may 
be succesfully reduced using the linearizing technique developed in 
this thesis. The compacting scheme for matrix equations effects gross 
savings for storage of large systems (Table I) and becomes even more 
attractive when the inverse of a matrix is either not desired or required. 
The Crank-Nicolson formulation coupled with Gauss-Seidel iteration as 
expressed in the computer code CRANKO has demonstrated itself to be a 
viable, extremely fast technique for the solution of differential 
equations. Further improvements to the CRANKO routine, specifically in 
regard to better control of the time step selection, may improve upon the 
attractiveness of this approach. 

Development is warranted in several areas. The computation of the 
Jacobian using iteration must be formulated in order to use the DVOGER 
subreutine with the compact system. DVOGER itself might usefully be 
reviewed to reduce its size requirements also. Investigation of the 
problems using the fine mesh size of 220 elements is continuing. This 


should yield more exact results and offer a better standard for comparison. 
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APPENDIX A 


TABLES 








oy ol EM 


NxN 
NxN 
Nxp 
NxN 
Nxp 


Nxp 


TABLE I 


IBM 360/67 CORE REQUIREMENTS 


in K Bytes 
NO. NODES OVERHEAD MATRIX STORAGE TOTAL 
38 100 225 6 PAS. 
38 100 O05) 135 
33 100 12 ey 
132 150 418 568 
Sys 150 38 188 
400 290 gba 402 t+ 


+calculated estimate 
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TABLE II 


COMPARISON OF METHODS 
38 NODAL PT. MODEL 


PROCESSING TIME (MIN) STEADY-STATE 
METHOD CORE (k CENTER UNIFORM OKEW ERROR (Z) 
DVOGER N° .01 330 > Bee Ga 7 
DVOGER N° .1 330 5.32 3.28 4.57 <.01 
DVOGER NxN exact 130 O22 6407 8.10 06 
DVOGER NxN approx 130 POR 3 Sloe S10h 7G 
CRANKO 110 62 20 22 3 
CRANKO DP 158 A, 23 25 5. 


Notes: 


ali 


Processing time is that time required by the program to reach one 
second in problem time. Steady state was obtained for each run. 


-O1 & .1 pertain to DVOGER error criterion EPS. All NxN trials 
used EPS = .1 CRANKO used iteration convergence tolerance of .0001. 


Double precision (DP) trials showed a difference in the fifth 
significant digit when compared to the SP counterparts. 


NxN approximate technique proceeds as rapidly as the we method 


until steady state is approached. This has been verified by 
experiments using various processing time cutoff points. 
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FIGURE 2 


41 
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COMPUTER PROGRAMMING CODES 


DISCUSSION OF PROCEDURES 
FED-2 LINEARIZED LISTING 
FED~2 COMPACT Nxp LISTING 
DATA PROCESSOR LISTING 


COMPACTING DEMONSTRATION LISTING WITH RESULTS 
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it. DISCUSSION 
A. The FED-2 Programs 

Both of the FED-2 programs consist of a small calling program which 
initiates the problem. The rest of the program is stored in the computer 
library. The advantages of a pre-compiled main program in time and 
convenience is clear, but the use of dummy dimension statements to 
pass data storage throughout the calling of the external subroutine 
DVOGER and its user-supplied subroutine appears as a trick not normally 
considered possible with the FORTRAN language. 

When the original i FED-2 program was written, every available 
device for reducing storage was attempted. Thus, where properties are 
normally summed over areas in the finite element method, the operations 
were carried out on the nodal points, since there are always fewer node 
points than elements. The new FED-2 programs continued these procedures 
not only to provide a valid comparison but also because this economizing 
technique is felt to be fundamentally sound even if unconventional. [Employ- 
ing this device, NFULEL is the number of node points associated with the 
core (fueled) region, whereas it would conventionally be the number of 
core (fueled) elements. 

For these programs, three interest points IPT1, IPT2, and IPT3, may 
be selected for detailed investigation, since the purview of information 
produced by the printed output of these programs is quite large, making 
detailed study difficult. The time history of the neutron flux (PSI) at 
each of these interest points is recorded throughout the problem in an 
array BATCH. If it is desired to use this collection, it can be written 
onto a computer storage by the inclusion of additional JCL in the GO step. 
In this manner, data was saved for analysis by the data processing program 


used in this thesis to prepare the graphs of Figures 4 - 21. 
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An options chart and a schedule for FED-2 input data deck listing 


are included in this section in addition to a listing of the programs. 


NCOUNT 
MTH 
ERRVAL 
A= 
ae 
EPSVAL 
1 
De 


mm O NM F 


OPTION SELECTIONS FOR FED-2 PROGRAMS 


PREMULTIPLICATION BY Vea 


PREMULTIPLICATION BY TRIAL 


DVOGER, Jacobian must be supplied 
DVOGER, computes its own Jacobian 
DVOGER, doesn't use Jacobian 


CRANKO routine specified, DVOGER not activated 


used internally by DVOGER, card input value 
has no effect 


used by CRANKO, it is the convergence criterion 
for satisfactory p values 


convergence criterion for DVOGER 


for CRANKO, sets the relative problem time 
difference required for a solution to be 
printed out 
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